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Abstract We propose a method for effectively upscaling incompressible viscous flow in 
large random polydispersed sphere packings: the emphasis of this method is on the determi- 
nation of the forces applied on the solid particles by the fluid. Pore bodies and their connec- 
tions are defined locally through a regular Delaunay triangulation of the packings. Viscous 
flow equations are upscaled at the pore level, and approximated with a finite volume numer- 
ical scheme. We compare numerical simulations of the proposed method to detailed finite 
element (FEM) simulations of the Stokes equations for assemblies of 8 to 200 spheres. A 
good agreement is found both in terms of forces exerted on the solid particles and effective 
permeability coefficients. 

Keywords viscous flow, granular material, solid fluid coupling, pore-network, finite 
volumes 



1 Introduction 

Understanding the mechanical behavior of fluid-saturated granular systems (e.g., soils, rocks, 
concretes, powders, and chemical reactors) hinges upon a detailed description of the stress 
exchange between fluid and solid phases. Fluid flowing through a matrix of solid spheres, 
exerts forces on the solid grains, displacing them from the trajectory they would have had if 
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List of Symbols 

a Non dimensional conductance factor 

£i Full domain of tlie two-pliase problem 

£ij Domain defined by tetrahedron i 

iiij union of tetrahdra (5,,,P,) and (Sij,Pj) 

r Part of Q occupied by the solid phase 

Fi Domain occupied by solid particle (' 

Part of Q occupied by the fluid phase (pore space) 

0j Part of £li occupied by the fluid phase (pore) 

0ij part of £2ij occupied by the fluid phase (throat) 

SiJ surface of the facet ij, separating tetrahedra i and j 

dx contour of domain X 

d^X pait of contour of X intersecting the fluid phase 

d'X part of contour of X intersecting (or in contact with) the solid phase 

Yij area of dOij in contact with spheres 

yfj area of the part of d0jj in contact with sphere k 

Sij the common facets of tetrahedra ^2, and Qj 

Ajj area of the fluid part S,y n of facet ij 

A^j area of the intersection Sij D -f|- of facet ij and sphere k 

Pj Voronoi dual (weighted center) of tetrahedra 

p' microscopic (pore-scale) fluid pressure 

Pi macroscopic fluid pressure in tetrahedra / 

u' microscopic fluid velocity 

u macroscopic fluid velocity 

V geometric contour velocity 

qjj flux through facet ij 

vf fluid volume contained in pore i 

Rjj hydraulic radius of throat ij 

R°ff effective radius of throat ij 

/J dynamic viscosity 

Ly length of throat ij 

f/ forces exerted by the fluid on the solid phase, x and y denote different terms in forces decomposition 

gij hydraulic conductance of facet (throat) ij 

Kij hydraulic conductivity of facet (throat) ij 



left to interact alone in a dry medium. In turn, the displacement of the grains creates a dy- 
namically changing domain for the fluid flow. The combination of these two effects results in 
nontrivial physically measurable effects, which are hard to capture by first-principles deriva- 
tions alone. Despite the important practical applications that are touched by fluid-saturated 
granular systems, no model currently exists that can efficiently reproduce the precise evolu- 
tion of large systems over different stress conditions in 3D problems. 

This work represents a first step in the direction of developing a fully-coupled, compu- 
tationally efficient model for the evolution of fluid-saturated granular materials under stress. 
In particular, we will focus our effort on the faithful approximation of the forces applied 
by the fluid on the solid grains, with the aim of incorporating these forces in the discrete 
element method (DEM) computations lCundall and Stracklll979l . 

While the movement of the solid grains can be efficiently modeled by DEM computa- 
tions, viscous fluid computations in dynamically changing domains present incredible com- 
putational challenges. At the microscopic (sub-pore) scale, fluid flow is governed by Stokes 
equations, which express fluid mass and momentum conservation at small Reynolds and 
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Stokes numbers: 



Vu = 



(1) 
(2) 



where u, and p are the microscopic fluid velocity and piezometric pressure, respectively, 
jU is the fluid dynamic viscosity, and 4" is a potential field (e.g., gravitational field). The 
piezometric pressure p is related to the absolute pressure p* via p = p* —p<P Equations 
([TJ-1I2I1 are augmented with a no-slip boundary condition for the fluid velocity, at the grain 
boundaries, u = 0, which is essentially responsible for the microscopic viscous energy losses 
(drag) that translate in a net loss of the macroscopic piezometric pressure, p, over the length 
of a porous column. 

At centimeter scale (and above), for relatively homogeneous porous materials, fluid flow 
is governed by Darcy law, which states the linear proportionality between the macroscopic 
pressure gradient Vp and the fluid flux q (discharge per unit area) : 



where k is the permeability. Darcy law emerges as a volume averaging of the viscous forces 
applied by the fluid flow at the surface of the grains over a sufficiently large number of pores. 

Recent advances in porous media imaging techniques give access to an unprecedented 
level of pore-space detail, down to the micrometer scale and below, where the Stokes equa- 
tions could be, in principle, solved numerically. The numerical solution of the Stokes equa- 
tions in spheres assemblies is, however, computationally expensive, and becomes prohibitive 
on commodity desktop computers for number of particles exceeding the hundreds. On the 
other hand, typical DEM simulations, can easily handle assemblies with number of particles 
ranging between 10'* and 10^ on the same com modity desktop comp uters, and compute their 
evolution over 10^ time-steps or more (see e.g. lScholtes et al.ll2009h . It is therefore obvious 
that a direct coupling of Stokes equations and DEM simulations becomes unfeasible for 
real-world problems, hence the need to develop an upscaled fluid model that minimizes the 
ratio between the degrees of freedom (DOFs) in the fluid and the number of solid particles. 

The past thirty years have seen a considerable effort in the exact upscaling of Equa- 
tions (|Tll-([2]i, i-S-, on how to obtain expressions for the permeability k from the knowledge 
of th e microscopic fl ow. Upscaling techniq ues such as the Volume Averaging with Clo- 
sure (IWhitakeilll999h . and homogenization dSanchez-Palencia and Zaouilll987h invariably 
require the solution of a (tensorial) Stokes problem, which has essentially the same algo- 
rithmic complexity of the "original" Stokes fluid flow problem. In this work we will not be 
concerned with this type of theoretical upscaling issues. Rather, we will focus our attention 
to the derivation of a physically based, simplified model of flow through porous materials 
that overcomes the numerical issues associated with the solution of the full Stokes problem. 

Different strategies have been adopted in the past for the solution of the fluid-solid cou- 
pling, which essentially differ in the modeling techniques adopted for the fluid part of the 
problem, as reviewed below. 

Microscale Stokes flow modeling 

The numerical solution of the Stokes equations is a computationally demanding task, espe- 
cially for complex three-dimensional pore geometries. Finite Element Methods (FEMs) are 
often used because of their flexibility in the definition of the numerical mesh. FEM meshes 




(3) 
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in three dimensions tend to be, however, very large and so do the computer memory foot- 
print and co mputational times need ed for the solution of the associated nonlinear system 
of equations. iGlowinski et aTl booih . for instance, reported computational costs of approxi- 
mately 0.1s per particle and per time-step in fluidization simulations on a SGI Origin 2000, 
employing an advanced fictitious domain method and a partially parallelized code: extrapo- 
lating these values to a number of particle typical of DEM problems, results in computational 
times in the order of 3-30 years. 

A numerical method that does not resort to the solution of large systems of nonlin- 
ear equations is the Lattice-Boltzmann (LB) scheme. While LB is generally faster than 
the FEM and has the possibility of being easily parallelized on multicore machines and 
CPUs (http://sailfish.us . edu.pl/ ), commonly implemented fixed size grids can re- 
sult in considerably larger computer memory occupancy in three dimensions. Only re- 
cently, grid refinement schemes have started to be implemented in open source LB codes 
(^tp : / /www . Ibmethod . org/palabos/ ). Realistic microscale flow simulations in com- 
plex pore geometries still require, however, access to large CPU clusters. 

Continuiim-discrele Darcy flow modeling 

In order to get acceptable c omputational costs, a number of authors cons idered c oarse-grid 
CFD methods (lNakasa"etal,.1999; Kawaguctii et alh998l : lMcNamara et al..2000uiKafui et alj 



l2002l : IZeghal and Shamvll2004l : lNiebling et al.ll2010h . In such methods, flow and solid-fluid 



interactions are defined with simplified semi-empirical models based on Darcy law. There is 
no direct coupling at the local scale: the forces acting on the individual particles are defined 
as a function of mesoscale-averaged fluid velocity obtained from porosity-based estimates 
of the permeability. Such strategies have been adopted for dilute particles suspensions as 
well as for dense materials at low strain rate (e.g., soils, rocks). Continuum-discrete cou- 
plings succeed in making coupled problems affordable in terms of CPU cost. The use of 
phenomenological laws for the estimation of the permeability, however, limits severely the 
predictive power of these models in uncalibrated parameter regions. Ultimately, such strate- 
gies do not render correctly the individual particle behavior, and cannot be reliably applied to 
problems such as strain localization, segregating phenomena, effects of local heterogeneities 
in porosity, and internal erosion by transport of fines. 

Pore-network modeling 

Pore-network modeling covers a third class of models, based on a simplified representation 
of porous media as a network of pores and throats. Pore-network models have been most 
con mionly developed to predict the permeability of materials from microstructural geome- 
try llBrvant et al.ll993l:lThompson and Foglei|l997l : iBakke and 0renlll997l : IPatzek and SilinI 



l200ll : lHilpert et al.ll2003l : lAbichou et al.ll2004ir but have also been extended to include mul 



tiphase flow effects (e.g., air bubbles, immi s cible two-phase flow ) llBrvant and Bluntll 19921 : 
lOren et al.iri998l : iBrvant and Johnsonll2003l : IPiri and Blunjlioosh . Crucial for their success 
is an adequate definition of how fluids are exchanged between pores in term of the local pore 
geometry. This aspect will be discussed further in section l2.2.2l 

Deformable pore-network modeling 

Pore-network models studies have mostly focused on flow in passive ligid solid frames. 
Little attention, however, has been devoted to the definition of forces applied to individual 
particles in the solid phase. The main aim of this work is to develop expressions for these 
forces, to be used in three-dimensional DEM simulations of solid particles immersed in a 
fluid flow. 
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Early ideas of a coupl ed pore-network flow and D EM can be found in the works of 
Hakuno and Tarumi ( Il995h and later by Bonilla (|2004) . These studies, however, were lim- 
ited to 2D models of discs assemblies, where pores were defined by closed loops of particles 
in contact. Since such pore geometry does not offer any free path for fluid exchanges, 2D 
problems implied some arbitrary definition of the local conductivity, assuming virtual chan- 
nels between adjacent voids. Adapting this approach to 3D spheres assemblies enables the 
definition of the local hydraulic conductivity using the actual geometry of the packing, as 
spheres packings always define an open network of connected voids. This in turns opens up 
the possibility of predicting both the macro-scale permeability and the forces acting on the 
individual particles rather than postulating it: this is the central theme of this work. 



2 Upscaled fluid flow model 

As described in the previous section, fine-scale Stokes flow numerical simulations present a 
prohibitive computational cost and an effective, upscaled fluid flow model is needed. In this 
section, we present a detailed derivation of the proposed upscaled fluid flow model, starting 
from the decomposition of the pore space in terms of pore bodies and local conductances. 
This decomposition results naturally in a finite volume description of the flow equations. The 
expressions for the forces on the individual particles are then derived. Finally, we describe 
the actual implementation of the model. 



2.1 Pore-scale volume decomposition 

Delaunay triangulations and their dual Voronoi graphs are widely used for st ructural studies 
of molecules, liquids, colloids, and granular materials (lAurenhammeijfl99ll) . More specifi- 
cally, it has been applied to dor nain decompositions in sphere packings, for the d efinition of 
microscale strains and stresses (ICalvettietal.ll9 97': Bagi 2006:'jeri er et al.l2010l). and pore- 
scale modeling of single-phase or m ulti-phase flow (.Brvant and Blunt! 19921 : Brvant and JohnsonI 



l2003l : lThompson and Fogle Jl997h . The most common type of Delaunay triangulation (here- 
after referred to as classical Delaunay) is defined for a set of isolated points. In this case the 
dual classical Voronoi graph of the triangulation defines polyhedra enclosing each point in 
the set, and facets of the polyhedra are parts of planes equidistant to two adjacent points. For 
mono-sized spheres, this construction gives a map of the void space. Classical Delaunay- 
Voronoi graphs, however, present a number of misfeatures when applied to spheres of dif- 
ferent sizes. Notably, branches and facets of the Voronoi graph can cross non-void regions, 
or divide regions in such a way that the splits are orthogonal to the actual void, as illustrated 
in fig.ffl 

iLuchnikov et al.|[l99^ proposed to define the split between spheres by means of curved 
surfaces equidistant from both spheres surfaces, thus defining Voronoi S cells. The use of 
Voronoi S graphs is, however, computationally expensive for large number of particles. 
Regular Delaunay triangulation (see fig.ffl, while overcoming the problems associated with 
classical Delaunay triangulations, provides a computationally cheap alternative to the use of 
Voronoi S cells. 

Regular Delaunay triangulation (lEdelsbmnner and Shahll993) generalizes classical De- 
launay triangulation to weighted points, where weights account for the radius of spheres. It 
can be shown that the dual Voronoi graph of the regular triangulation is entirely contained in 
voids between spheres, as opposed to the dual of classical Delaunay, as seen on fig.ffl (2D) 
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Classical Delauray Regular Delaunay 




triangulation triangulation 

Voronoi graph Voronoi graph 



Fig. 1 Comparison of 2D triangulations and dual Voronoi graplis obtained from classical Delaunay and reg- 
ular' Delaunay: classical Voronoi graph has branches inside discs, while regular Voronoi gives all branches in 
the voids space. Observe that if R tends to infinity, {P[,Pi,P2,P2} tends to aligned points. 




Fig. 2 Adjacent tetrahedra in the regular Delaunay triangulation and dual Voronoi network, in two dimensions 
{a,b) and three dimensions (c). 

and fig.|2](2D/3D). Edges and facets of the regular Voronoi graph are lines and planes, thus 
enabling fast computation of geometrical quantities. In section l2".2.2l we will use a combina- 
tion of regular Delaunay facets and regular Voronoi vertices to decompose the pore volume. 

In what follows, Q denotes a domain occupied by a porous material, in which F and 
are the domains occupied respectively by the solid and the fluid: Q = FU0,Fn0 = 0(0 
is also called pore space). We denote by Nc the number of tetrahedral cells in the regular 
Delaunay triangulation of the sphere packing, and Q, the domain defined by tetrahedron ; : 
£2 = u|^ji2/. Similarly, Ns is the number of spheres, and /] the domain occupied by sphere 
i, so that F = ufj^^Fi. 



2.2 Fluxes 
2.2.7 Continuity 

We define 0, as the portion of the tetrahedral element X2, not occupied by spheres. The 
volume vf of 0, is filled by fluid (fig.[3]l, since we are considering saturated porous media. 
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Fig. 3 Tetrahedral element of the finite volume decomposition. 

The continuity equation for an incompressible fluid (|2}, cast into its surface integral form 
and using the divergence theorem, gives a relation between the time derivative of V- and the 
fluid velocity: 

v[ = j (v-u)-ndi, (4) 

50, 

where 50, is the contour of 0,, n the outward pointing unit vector normal to (90, u the fluid 
velocity, and v the velocity of the contour. One part of the contour, noted (9'0, corresponds 
to a solid-fluid interface. At any point in (9*0, we have (v — u) ■ n = 0, so that the integra- 
tion domain above can be restricted to 0, the fluid part of the contour. By introducing 
sjjii G {71572,73,74} the intersections of triangular surfaces Sjj with the fluid domain (fig. 

[3]l, so that 0/ = ^f^ji^ij, we can define four intergrals describing fluid fluxes qij from 
tetrahedron / to adjacent tetrahedra ji to 74. Finally: 

v/= t /(v-u)-nd.= - (5) 

In deformable sphere packings, Vj can be computed on the basis of particles motion (i.e. 
velocities of the vertices of the tetrahedra), thus linking fluid fluxes with the deformation 

of the solid skeleton. The derivation of V,- as function of particles velocities is not detailed 
here for brevity, because in what follows particles are assumed to be fixed: Vf is constant 
and the continuity equation simplifies to 

E'?'7 = 0- (6) 

j=ji 



2.2.2 Local conductance 

Both Stokes ([T]l and Darcy ([3]l equations imply a linear relation between pressure gradients 
and fluxes. Here, we introduce a inter-pore gradient, defined as the ratio between pressure 
difference pi — pj between two connected tetraheral cells, and a relevant length L,j - to 
be defined below. Being linear, the relation between cjij and the inter-pore gradient can be 
expressed using the local conductance factor gij of facet ij: 



(7) 
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Regular Triangulation 

Regular Voronoi Tesselation 

I I !!,; 

an, 
— fiij 
111 e„ 



Fig. 4 Construction of subdomain ©ij, defining the volume of the throat assigned to a facet for the definition 
of hydraulic radius R*.; in 3D (left) and 2D (right). 



Considerable efforts have been devoted in the past literature to the definition of g,y in pore- 
network models, with attempts to generalize Poiseuilles's law to pores of different shapes, 
eventually mapping the microstructure o f some real materials for permeability predictions 
jPiri and Blun3f2005l : Iffilpert et al.ll2003l) . By analogy with the Hagen-Poiseuille relation, 
gij may be defined by introducing the hydraulic radius of the pore throat R'Ij, and its cross- 
sectional area Ajj (the definitions of Rjj and A,y are discussed below), by means of a non- 
dimensional conductance factor a reflecting the throat's shape (Hagen-Poiseuille for circu- 
lar cross-sections of radius IR'fj is recovered with a = 1/2): 



8U ■ 



(8) 



Lij, R'Ij, andAij are geometrical quantities describing the throats geometry. Even though 
these variables are found in most of the papers cited above, there is no general agree- 
ment on their definition for pore-network modeling of arbitrary microstractural geometry. 
A detailed analysis of the effect of non-circular cross-sections with variable cons trictio n 
along the flow path can be f ound i n Patzek a nd Silin 2001 and JVlortenseri et . aO 1 2005h . 
while [Thompson and Fogled ( ll997hlBakke and 0renl ( Il997lb . and iBrvant et air ( ll993l) fo- 
cused specifically on triangulated sphere packings. Some of these models, however, do not 
clearly define a partition of the void space, and as such they need ad-hoc corrections for L,j 
to ensure that the s ame pore volumes are not accounted for multiple times in different cells 
tervant et all 19931) . We observe here that in classical Delaunay-Voronoi graphs of spheres 
with polydispersed radii , circumcenters (and barycenters) may lie inside the solid phase 
tervant and BluntI ( Il992h ). and definitions of L,j based on distances between circumcenters 
become problematic. 

For the present model, we take advantage of the regular Voronoi graph structure, whose 
edges and vertices are always contained in the pore space. A partition of the full domain is 
defined as the set of sub-domains Qjj, with Qjj the union of two tetrahedra constructed from 
facet Sij and Voronoi vertices P, and Pj : Qij = {Sij,Pi) U {Sij,Pj). Keeping only the part of 
Qij that intersect the pore space, we define a throat connecting two pores, noted 0,j (see 
figs. |4] and [Sj. We define the hydraulic radius R'Ij in eq.[8]as the ratio between the throat's 
volume and solid-fluid interfaces area. Denoting by the volume of 0^, and yij the area 
of d''@jj (the part of the contour in contact with spheres), the hydraulic radius R'-: reads: 



(9) 
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Note that, with the partition of the voids space proposed here, the definition of effective radii 
accounts for the exact total solid surface and void volume in the packing. Noting ^ the set 
of facets, we have = IJijg^ ^nd d© = U/ye,^ d^Qij. 

We define the length of the throat as the distance between Voronoi vertices: L,j = 1| , 
and the throat's cross-sectional area as that of the fluid domain Sjj introduced in the previous 
section (fig.[3]l. In this model, the factor a is uniquely assigned to all throats of a packing 
and can be considered a calibration factor. A value of a = 1/2 has been choosen as an initial 
guess, by analogy with the Hagen-Poiseuille law. The numerical simulations presented in 
section [3] indicate that setting alpha=l/2 is indeed appropriate for the spheres geometries 
considered in this investigation. We observe that the values of a obtai ned analytically v ary 
in a moderately wide range of values for compact throats sections (iPatzek and Silinl200lh . It 
is also worth noting that, the factor a being the same for all throats, modifying this value will 
not modify the pressure distribution in the boundary value problems presented below (mixed 
Neumann-Dirichlet), and consequently it would not affect the values of forces applied on 
the particles. 

Another possible definition of gij has been proposed by iBrvant and Bhin^ (Il992h for 
dense packings of mono-sized spheres. These authors defined gij as a function of an ejfective 
radius R'jj^ (ER), accounting thus only for the effects of the narrowest cross-section (surface 
Sjj in fig.lHl and disregarding the influence of the full volume of the flow path. The effective 
radius is defined ' as R^j^ = {/J' + r'"j^'^''-)/4, where rJJ' is the radius of the disk of same 

surface as sfj, and r™"''' is the radius of the circle inscribed in the three spheres intersecting 
the facet. It represents the hydraulic radius of the circular tube that would have the same 
conductance as the throat: 



271 R'V^ A'/Zr'//^ 



(10) 



M 2m ' 

where AJJ-^ = 7t{2R'jj-^)^ is the cross-sectional area of the equivalent tube. The second ex- 
pression is given for the purpose of comparison with eq. |8] The values of R'^j and R^j^ in 
random dense packings of polydispersed spheres, and the results they give in terms of fluxes 
and forces are compared in section |3] 



2.3 Forces on solid particles 

The total force generated on particle k by the fluid includes the effects of absolute pres- 
sure p* and viscous stress T, which divergence is the right-hand side of ([T): 

F*= j {p*n + Tn) ds. (11) 

We recall that the piezometric pressure governing the flow problem was defined as p = 
p* — p4>(x). By introducing p in the previous equation, one can define F*^ as the sum of 
three terms, noted F* •^ F'^ *^ and F'' *^: 

F'^= J p4>(x)ndi+ J pnds+ J rnd.? = F" -FF^-* -FF"-*. (12) 

an 3n an 

' For the convenience of comparisons, and consistently witli tfie expression of conductance we are using, 
we introduce the effective radius R'' f f as the half of its definition in the original paper from Bryant and Blunt. 
The final value of the conductance is preserved. 
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Voronoi graph 
Delaunay graph 
Cell dual 



Fig. 5 Volume decomposition for viscous drag force definition : (a) pressure distribution on 30,j (in 2D for 
clarity), (b) definition of facet-spheres intersections in 3D. 



F* * is the so-called buoyancy force, and can be computed independently. In the case of 
gravitational body forces, it will give the Archimede's force proportional to the volume 
of F''"* and F''* are terms resulting from viscous flow: F'' *^ resulting from losses of 
piezometric pressure, and F''* resulting from viscous shear stress. Both F'' * and F'' * will 
be derived separately, at the scale of the domains Qjj that were introduced in the previous 
section, in the approximation of piecewise constant pressure. 



2.3.1 Integration of pressure 

The force generated by p on particle k in domain Qjj is the sum of two terms implying 
pressures pi and pj (see fig.O, 

F^}'*= J pinds+ J pjnds (13) 

In three dimensions, computing such integral on spherical triangles could be computation- 
ally expensive. It is more convenient to project the pressure on the conjugate planar parts 
(angular sectors) of the closed domain /| OdQij, which trace in the plane of fig. |5] corre- 
sponds to segments OO' and OO". Note that when iterating over all domains Qij adjacent 
to one particle, the integral on sectors of type OO" will appear twice with opposite normals 
and cancel out each other. Finally, the contribution to pressure force on particle k in domain 
Qij is simply proportional to the area of sector Sij nil (trace 00'). If n,y is the unit 
vector pointing from P, to Pj, the force reads: 

Fff=4{pi-pj)nij (14) 

2.3.2 Integration of viscous stresses 



To integrate the viscous stresses, we first define the total viscous force FJ^ applied on the 
solid phase in Qjj. Since Qij intersects three spheres, F,y will have to be later splited into 
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three terms. F,y is defined as 



F;}= J tnds. (15) 



An expression of FJ^ is obtained by integrating the momentum conservation equation 
lUll) in 0ij. The integral is cast in the form of a surface integral on contour d@ij using the 
divergence theorem : 

j (pn + Tn) di = 0. 

50,, 

This integral can be decomposed as in eq.[T6l where the terms correspond respectively to FJ^ , 
to the sum of viscous stress on the fluid part @ij of the contour, and to the sum of pressure. 
By neglecting the second term (thus assuming that pressure gradients are equilibrated mainly 
by viscous stress on the solid phase), the expression of F,y takes the form of equation [17] 
where the viscous force is simply proportional to the product of the throat's cross-sectional 
area Ajj and the pressure jump pj — p,. 

Tnd.?-!- 1 Tnd.v+ j pnAs = Q. (16) 

Y]j ^ - I pn d^ = a{. {pj - Pi) mj (17) 

In order to define the viscous forces applied on each of the three spheres intersecting 
Qij, it is assumed that the force on sphere k is proportional to the surface of that sphere 
contained in the subdomain. If yfj denotes the area of the curved surface OQij, the force 
on k then reads: 

r:f = F:,-fi^ (18) 

Finally, the total force on one particle is obtained by summing viscous and pressure forces 
from all incident facets with the buoyancy force: 

L {F;;*+Ff/}+F^'* (19) 

('7) incident 



2.4 Implementation 
2.4.1 Network definition 

The network model has been implemented in C ++, and it is freely a vailable as an op- 
tional p ackage of the open-sou rce software Yade dSmilauer et alj|2010l) . The C++ library 
CGAL teoissonnat et alj|2002h is used for the triangulation procedure. This library ensures 
exact pre dicates and c onstructions, and is one of the fastest computational geometry codes 
available jLiu and Snoevink .20051). It is worth noting that regular triangulation involves only 
squared distances comparisons, thus avoiding time consuming square roots and divisions. 
The computation of forces (eq. [T9t also requires only simple vector multiplications. The 
only non-trivial operation is the computation of solid angles needed to define the volumes 
and surfaces associated to sphere-tetrahedron intersections, in eq|9l Th is cost can be reduced 
significantly, however, using the algorithm of Oosterom and Strackee ( Il983h . 
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2.4.2 Flux problem solution 

Combining Equations |6] and [7] gives for element i: 

E J^{Pi-Pj)= I Kij{pi-pj) = v[ (20) 

j=l,..,4^'V j=l,..,4 

If the triangulation defines A^^. tetrahedra and if the number of elements with imposed 
pressure is A',,, the pressure field is obtained after solving a linear system of size Nc — Np. 
The matrix of the system is sparse and symmetric {Kij = Kji), and positive defined. An over- 
relaxed Gauss-Siedel algorithm has been used to solve this linear system in the simulations 
presented below. This algorithm has been chosen for its simplicity, though other strategies 
could be preferred in the future for better performance and for complex couplings (including 
particles motion). 



3 Comparison with small-scale FEM simulations 

3.1 Numerical setup 

In this section, we compare the results of the pore-scale finite volume (FV) formulation 
described in section |2] with FEM Stokes flow simulations of Equ. ([Tli-([2]i, for dense sphere 
packings subjected to an imposed pressure gradient. We consider sphere packings contained 
in a cube of size Iq (all results below will be normalized with respect to this reference length). 
Boundaries are accounted for in the triangulation process in the form of spheres with near- 
infinity radii (10* x Iq). Hence, planes are not introducing special cases to handle in the 
algorithms, and equations presented in previous section apply for boundaries as for any 
other sphere. 

The smallest test problem we consider is a regular cubic packing of 8 spheres. Packings 
of nine spheres are obtained by placing an additional sphere (with variable size d) in the 
center of the cubic packing (fig. |6}. Larger assemblies of 25, 54, and 200 spheres (fig. [7) 
are ra ndom packings that were generated using the DEM software Yade jSmilauer et al.l 
I2OIOI) by simul ating t he growth of spheres after random positioning, using the algorithm of 
Chareyre et al. Radii are generated randomly with respect to uniform distributions 

between radii ?-„„„ and r,„ax- In the 200-spheres assembly, sizes are narrow graded, with 
rmin/rmax = 0.9 (this Small dispersion of sizes prevents the formation of crystal-like patems). 
Sizes in the 25-spheres and 54-spheres packings, are more dispersed, with r„,i„/r„,ax = 0.5. 
After dense and stable packings are obtained, the spheres's positions are fixed. 

For each packing geometry, the flow boundary conditions are the ones described on 
fig. |6l pressure is imposed on the top and bottom boundaries, while no-flux conditions are 
imposed on the lateral boundaries. 

Both "no-slip" and "symmetry" conditions are considered on the lateral boundaries. In 
the FV modeling, they are reflected in the definition of conductivity and forces in the throats 
©ij in contact with boundaries as follows: 

- For the "no-slip" condition, the surfaces of infinite spheres in contact with 0,y are ac- 
counted for in Eqs.|9]and[T8] as for any other sphere. 

- For the "symmetry" condition, these surfaces are ignored : )fj = 0, when k is the indice 
of a bounding sphere. 



p = 1 




Fig. 6 Boundary conditions of flow simulations (a), FEM mesh (b) and FV mesh (c) of a 9-spheres packing. 
The FV packings are plotted with Voronoi graph, where pressure values are defined at each point. 




Table [T] gives a comparison of problems sizes in terms of degrees of freedom (DOF) of 
the pressure field, and CPU times for solving. The data in this table is only an indication of 
how much is gained from the pore-scale formulation of Stokes flow: impacts of coarsening 
the FEM mesh, or benefits of other algorithms in the FV model (e.g. conjugated gradient) 
have not been investigated yet. 
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Table 1 Compared DOF's and CPU time, (a) includes only solving Stokes flow, excludes preprocessing (e.g. 
mesh generation) and postprocessing (forces on spheres) done via a graphical user interface; (b) includes 
packing triangulation, solving, and computation forces on particles; (c) out-of-memory; no result. 



Nb. of 


FEM dofs 


FV dofs 


FEM 


time'"' FV time'*) 


spheres 






(s) 


(s) 


9 


1.7e5 


45 


300 


0.00022 


200 


1.2e6 


1093 


5400 


0.0046 


2e3 


(c) 


11.6e3 


(c) 


0.091 


2e4 


(c) 


11.3e4 




2.21 




Fig. 8 Isovalues of pressure in the 200 spheres packing on the plane .v = 0.5/o, obtained with FEM (a) and 
FV (b). 



3.2 Pressure Field 

The pressure fields from FEM and FV are compared on figure |8] Local conductances are 
defined as in eqs. l[8]|9ll. The isovalues in the FV result are, by model definition (element- 
wise constant pressure), coincident with facets. Thus, the curvature of some isolines, as 
observed in FEM results, cannot be reproduced. The curvature of the isolines in the FEM 
solution is, however, usually very small so that most of them are reasonably approximated 
by straight lines. Furthermore, most isolines are found near the necks of the flow paths, thus 
justifying a posteriori the approximation of element-wise constant pressure used in the FV 
scheme. Overall, the two fields compare well. The deviation from the horizontal isolines 
that would be obtained in an homogeneous continuum with same boundary conditions is 
relatively well reproduced by the FV model. 



3.3 Fluxes 

Two fluxes can actually be obtained from simulations: influx Q,, and outflux Qo at the re- 
spective boundaries (in the FV model, they are computed by summing the qij of eq.[8]for 
elements where pressure is imposed). It was found that the difference between g, and Qo 
is always negligible (less than 10^^ x |2,|). This is consistent with the incompressibility 
assumption and indicates a good convergence of the numerical solvers, be it in FEM (direct 
sparse solver Pardiso) or FV (Gauss Siedel method). In the following, we define Q = Qj. 
Fluxes are expressed in the form of normalized permeability defined as K = ^ Qho/{APS^), 



15 




Number of spheres 



Fig. 9 Predicted permeability in FEM and FV versus size of packing, for no-slip (NS) and symmetric (sym) 
boundary conditions. The FV results include conductances defined using the hydraulic radius (HR - eq. |9) 
and effective radius (ER - eq. llOt . 



with S = the packing's cross-sectional area. Note that the introduction of permeabihty is 
only for the convenience of results representation, and is not implying anything with respect 
to the representative volume element (REV) concept. The topic of this paper is the compar- 
ison of small-scale modelings of fluxes and forces for some spheres packings: derivation of 
higher-scale properties of an equivalent continuum is beside the scope of the present inves- 
tigation. 

Permeabilities obtained with FEM and FV for the different packings are compared on 
fig. 121 Both no-slip and symmetry conditions are considered for lateral boundaries. Compar- 
ison between the definitions of conductivity using the "hydraulic" radius (HR) of eq.|9]and 
Bryant's "effective" radius (ER) of eq.[TO]are provided. This later comparison only covers 
no-slip conditions, as the definition of effective radii in the symmetry case was ambiguous. 

Symmetric boundary conditions give an average ratio of 1.01 between Khr and Kpem, 
with a maximum deviation of +13% for the 25-sphere packing. With no-slip conditions, 
the ratio are 0.78 for Khr/Kpem (max. deviation —40%) and 0.39 for Khr/Kfem (max. 
deviation —63%); Khr is giving the best estimate of Kpem in all cases. The fact that the FEM 
results are better reproduced for symmetric conditions suggests that the fluxes along planes 
are underestimated. We consider, however, that the predictions of fluxes from hydraulic 
radius can be considered satisfactory overall. 

The evolution of permeability in 9-spheres packings as a function of size d of the inner 
sphere's size (fig.|6]l is plotted on fig. [TO] The evolution of K with d/D is correctly reflected 
in HR-based results. Again, FV results match FEM better for symmetric boundaries. ER- 
based simulations are underestimating the fluxes by an average factor Ker/Kfem = 0.43. 

It can be concluded that the initial value a = 1/2 of the conductance factor entering eq. 
[8] tends to underestimate fluxes in average. Although systematic comparisons could enable 
a better calibration of a, deviation from the FEM results is an inherent defect due to the 
approximations of the pore-scale description adopted here. Considering the small number 
of samples we tested, a has not been further adjusted to closely match FEM results. 

To investigate further the differences found between Khr and Ker, the values of and 
R^f f obtained for all facets of the 200 spheres packing triangulation have been plotted (fig 
fTTT i. Clearly, the two values are strongly correlated, with the exception of very few points 
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FEM no slip 
FEM symmetry - 
FV/HR - no slip 
FV/HR - symmetry 
^ FV/ER - no slip - 



0.2 0.3 



Fig. 10 Permeability of 9 spheres packing versus size of the the central particle. 
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Fig. 11 Effective radius versus hydrauUc radius for all facets in the 200 grains packing (a), with exemple of 
a facet giving S^J' > S*^ . 



that will be commented below, but they deviate from the x = y line. The average value of 
and R'ff is 1.5. Considering this distribution, it is obvious that the differences found 
between Khr and Ker are reflecting the fact that 7?'' > R'^^^ in average at the local scale. 

A small number of points fall out the well correlated cloud, for which R'^f f >/?''. A close 
inspection of the outliers revealed that these correspond to comer cases similar to the one 
represented on fig.[IT] where a flat triangle results in the inscribed circle overlapping outside 
its original facet. In such a situation, a flatter triangle will result in a smaller hydraulic radius 
(smaller but a higher effective radius. The inscribed circle will eventually overlap other 
inscribed circles or spheres of the packing, thus loosing physical consistency: it was not 
clear to us how such special cases should be handled. The number of affected facets is, 
however, small and do not significantly affect the results. These results lead us to conclude 
that he hydraulic radius is overall a more robust parameter than the effective radius. 
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Table 2 Nomialized forces in 9 spheres packings (d/D = 0.72). All forces ai'e divided by SAP, so that the 
sum of all forces should be exactly 1 . 





FEM 




FV 






Forces on one big sphere (y-component) 










viscous force 


1.51 X 10- 


-2 


1.06 > 


c 10- 


-2 


pressure force 


9.09 X 10" 


-2 


9.64 X 10" 


-2 


total 


1.06 X 10- 


-1 


1.07 > 


c 10- 


-1 


Forces on the small sphere (y-component) 










viscous force 


9.38X 10- 


-3 


5.83 > 


c 10- 


-3 


pressure force 


5.17 X 10- 


-2 


5.72 > 


; 10- 


-2 


total 


6.04 X 10- 


-2 


6.30 > 


; 10- 


-2 


Forces on boundary x=0 












viscous force (y-component) 


2.27 X 10- 


-2 


2.03 > 


c 10- 


-2 


pressure force (x-component) 


4.98 X 10- 


1 


5.00 > 


c 10- 


1 


Total force on the solid phase 












y-component 


1-5.56X 10-* 


1 + 1.18X 10-'' 



3.4 Forces 

Forces on particles and boundaries, as defined in section 12.31 have been obtained for all 
packings using the FV scheme. In all cases, the sum of pressure and viscous forces on 
spheres and boundaries is, as expected, close to AP x S (±10^^%), where S is the packing's 
cross-sectional area. Comparison for the 8-spheres packing is trivial: Fi = S x AP/% on each 
sphere for both FEM and FV, modulo numerical errors. A detailed comparison is presented 
here only for the 9-spheres packing, for which contour integral of fluid stress in the FEM 
models (the FEM results in table|2]have been collected. 

Table [2] gives the total fluid forces per element of the solid phase, and details the vis- 
cous and pressure contributions. The viscous force on lateral boundaries, total force on big 
spheres, and total force on the center sphere, are all approximated with an error smaller that 
10%. Surprisingly, the prediction is less good for individual viscous and pressure contribu- 
tions than for the total force. The viscous part is overestimated, while the pressure effect is 
underestimated. The two errors compensation results in a smaller error on the total force. 

The evolution of force terms and error for the center sphere is given on fig. [TJ] This 
result shows an increasing error with decreasing d/D. For d/D > 0.1, the FV scheme un- 
derestimates the force on the the center sphere by a factor 2. This suggest a limitation of the 
current formulation in the limit of (1) high contrasts in particles size, and/or (2) small par- 
ticles "floating" in the voids between big ones. Here again, more investigations are needed 
in order to determine what situation exactly is generating this deviation from the FEM so- 
lution. As long as d/D > 0.2, though the predicted force is very well predicted, with error 
less than 10%. It is observed that similar forces are obtained with hydraulic radius (HR) and 
effective radius (ER). 

4 Conclusion 

The pore-space of a dense spheres packing can be efficiently represented by means of a reg- 
ular Delaunay triangulation. In this way, the total number of DOFs and CPU time required 
for calculating flow and forces acting on the particles are reduced drastically with respect to 
small scale Stokes flow FEM calculations. 
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FEM (total) 
FV/HR (total) — a— 

FV/ER (total ,g - 0.5 




d/D 



Fig. 12 Total force and viscous force applied on a particle of size d, placed in the center of a cubic packing of 
8 particles of size D, in FEM simulations, FV simulations based on hydraulic radius (HR), and FV simulations 
based on effective radius (ER). 



Expressions of local conductivities and forces induced on the particles have been de- 
rived. Bounding planes are straightforwardly taken into account in the method, and their 
effect in terms of permeability and viscous forces is reflected in the model for both no-slip 
and symmetry conditions. The method could be extended to account for non non-spherical 
particle shapes, as approximated by the intersection of a finite number spheres of various 
diameter. 

The definition of local conductivity of pore throats is based on a definition of the local 
hydraulic radius R'', as in the Kozeny-Carman (KC) derivation of permeability in porous 
material. The difference with KC-like derivations lies in the fact that connectivity and tor- 
tuosity of the network do not have to be included in the derivation, since they arise directly 
from the tetrahedral domain decomposition. The definition of the hydraulic radius R'' is 
simple and involves almost only real and 3D vectors multiplications for the definition of 
surfaces and volumes, which favors numerical performance. 

Taking the FEM results as reference, the permeabilities obtained with the pore-scale 
modeling are in most cases within the range of ±20%, usually below the reference value. 
Considering only symmetric boundary conditions, the permeabilites are better predicted (in 
the range ±10%), which indicates that fluxes along planar no-slip boundaries are under- 
estimated. We consider that such estimates are acceptable overall, keeping in mind that 
permeability prediction is of secondary importance, in comparison with forces, for the hy- 
dromechanical couplings. 

The effective radius proposed by Bryant et al. (Il993h for mono-sized spheres also gives 
acceptable predictions of permeabilities for polydispersed sphere packings, but it is com- 
putationally more expensive and underestimates fluid fluxes more than the hydraulic radius 
proposed in this investigation. It is important to note that the pore network model of Bryant 
et al. includes a reduction of , accounting for overlapping cylindrical throats. This correc- 
tion has not been implemented because it needs relatively complex computations of cylin- 
ders intersections. By maximizing micro-gradients, this correction could have at least partly 
balanced the underestimation of permeabilities in our results. Hence, our findings cannot be 
considered in disagreement with Bryant et al. model in itself. They only indicate that the 
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effective radius should not be used in connection witii the regular Delaunay partitioning we 
are developing. 

The pressure field compares well with the FEM results. In all cases, the sum of forces 
applied on the solid phase in a unit cube under unit gradient is close to one, showing the cor- 
rect implementation and validity of the model at the mesoscale. Forces applied by the fluid 
on individual particles, and viscous forces on no-slip boundaries, are found to be correct, 
with errors typically less than 10%: the error, though, grows larger when the size ratio d/D 
of the particles is less than 0.2. 

This modeling of fluid forces gives a sound basis for fluid-particles systems. Current 
works focus on the fully coupled problems in deformable packings with deformation-induced 
pore pressure, and on the optimization of the flow computation in this context. 
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